import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import seaborn as sns

def perform_monte_carlo_clustering(file_path, n_simulations=1000, sigma=0.0042, k_clusters=5):
    """
    多次扰动模拟 Y 染色体浓度，进行 KMeans 聚类，并统计 BMI 分组和聚类中心
    """
    sns.set(style="whitegrid")

    # 1. Load data
    df = pd.read_excel(file_path)
    bmi_col = '基准BMI'
    y_col = 't90'  # 假设 t90 是 Y 染色体浓度

    # 2. 保存每次模拟的结果
    cluster_centers_all = []
    bmi_ranges_all = []

    for sim in range(n_simulations):
        # 2a. 添加扰动
        noisy_y = df[y_col] + np.random.normal(loc=0, scale=sigma, size=len(df))

        # 2b. 聚类
        X = noisy_y.values.reshape(-1, 1)
        kmeans = KMeans(n_clusters=k_clusters, init='k-means++', random_state=42, n_init=10)
        clusters = kmeans.fit_predict(X)
        df['cluster_sim'] = clusters

        # 2c. 保存聚类中心
        cluster_centers_all.append(sorted(kmeans.cluster_centers_.flatten()))

        # 2d. 保存 BMI 范围
        bmi_ranges = df.groupby('cluster_sim')[bmi_col].agg(['min', 'max']).reset_index()
        bmi_ranges_all.append(bmi_ranges[['min', 'max']].values)

    # 3. 转为 NumPy array 方便统计
    cluster_centers_all = np.array(cluster_centers_all)  # shape: (n_simulations, k_clusters)
    bmi_ranges_all = np.array(bmi_ranges_all)  # shape: (n_simulations, k_clusters, 2)

    # 4. 统计每个聚类中心的均值和标准差
    cluster_centers_mean = cluster_centers_all.mean(axis=0)
    cluster_centers_std = cluster_centers_all.std(axis=0)

    print("每个聚类中心均值与标准差:")
    for i, (mean, std) in enumerate(zip(cluster_centers_mean, cluster_centers_std)):
        print(f"聚类 {i}: 均值={mean:.4f}, 标准差={std:.4f}")

    # 5. 统计 BMI 范围
    bmi_min_mean = bmi_ranges_all[:, :, 0].mean(axis=0)
    bmi_min_std = bmi_ranges_all[:, :, 0].std(axis=0)
    bmi_max_mean = bmi_ranges_all[:, :, 1].mean(axis=0)
    bmi_max_std = bmi_ranges_all[:, :, 1].std(axis=0)

    print("\n每个分组 BMI 范围均值 ± 标准差:")
    for i in range(k_clusters):
        print(f"分组 {i}: BMI_min = {bmi_min_mean[i]:.2f} ± {bmi_min_std[i]:.2f}, "
              f"BMI_max = {bmi_max_mean[i]:.2f} ± {bmi_max_std[i]:.2f}")

    # 6. 可视化聚类中心的分布
    plt.figure(figsize=(10, 6))
    for i in range(k_clusters):
        sns.kdeplot(cluster_centers_all[:, i], label=f'聚类 {i}')
    plt.xlabel('聚类中心 (t90)')
    plt.ylabel('密度')
    plt.title('多次扰动模拟后的聚类中心分布')
    plt.legend()
    plt.savefig('cluster_centers_distribution.png')
    plt.close()

if __name__ == '__main__':
    perform_monte_carlo_clustering('data_after_q2_2.xlsx', n_simulations=1000, sigma=0.0042, k_clusters=5)
